import numpy as np
import matplotlib.pyplot as plt
import pickle
from astropy.io import fits
import os
import sys
import copy

from fit_HDO_H218O_template_mcmc_final import *

kb=1.36e-16
hp=6.626e-27
c=2.99792458e10
d=400.0
pc=3.09e18
jy=1.0e-23
pi=3.14159
B_0=272912.6e6
A_ul_225=1.319e-05
A_ul_241=1.187e-05
A_ul_203=4.813e-06

E_u_225=167.56087 
E_u_241=95.22422  
E_u_203=203.7

g_u_225=7.0
g_u_241=5.0
g_u_203=7.0

nu_225=225.896720e9
nu_241=241.561550e9
nu_203=203.40752000e9

dv_241=0.4
dv_225=0.4
dv_203=0.4

def jy_beam_to_k(bmaj=None,bmin=None,freq=None):
    import math
    omega=(pi*bmaj*bmin)*(pi/180.0)**2
    wave=c/(freq)
    jy2k=wave**2*jy/(2.0*kb*omega)*4.0*math.log(2.0)
    return jy2k

def jy_to_k_aperture(bmaj=None,bmin=None,freq=None):
    import math
    omega=(pi*bmaj*bmin)*(pi/180.0)**2
    wave=c/(freq)
    jy2k=wave**2*jy/(2.0*kb*omega)
    return jy2k

def jy2k(freq):
   wave=c/freq
   return jy*wave**2/(2.0*kb)

def Qrot_HDO(T_ex):
   Q=np.array([1.334,2.711,6.952,18.850,52.300,95.565,146.859])
   T=np.array([9.375,18.75,37.5,75.0,150.0,225.0,300.0])
   Q_interp=np.interp(T_ex,T,Q)
   return Q_interp

def Qrot_H218O(T_ex):
   Q=np.array([1.26,03.054,8.648,23.361,64.210,117.004,179.639])
   T=np.array([9.375,18.75,37.5,75.0,150.0,225.0,300.0])
   Q_interp=np.interp(T_ex,T,Q)
   return Q_interp

#provides total number of molecules, not a column density, needs erg cm^-2 s^-1 input
def N_u_Jy(intensity,nu_u,A_ul):
  gamma=4.0*3.14159*(d*pc)**2/(hp*nu_u*A_ul)
  N_u=intensity*gamma
  return N_u

#provides molecules per unit area, needs K*km/s input
def N_u_K(intensity,nu_u,A_ul):
   gamma=8.0*3.14159*kb*nu_u**2/(hp*c**3*A_ul)
   N_u=intensity*gamma
   return N_u

import sys

if len(sys.argv) != 2:  # the program name and the two arguments
  # stop the program and print an error message
  sys.exit("python diskflux.py measurement_type (H218O, methanol, mixed, kepmask, kepmask-conv, kepamsk-ul, or stacked-methanol)")

measurement_type = str(sys.argv[1])
if measurement_type != 'H218O' and measurement_type != 'methanol' and measurement_type != 'mixed' and measurement_type != 'kepmask' and measurement_type != 'kepmask-5-to-7' and measurement_type != 'kepmask-5-to-7-scale' and measurement_type != 'kepmask-ul' and measurement_type != 'mom0'  and measurement_type != 'mom0-scale' and measurement_type != 'stacked-methanol':
   sys.exit('invalid measurement type selected please choose (H218O, methanol, mixed, kepmask, kepmask-conv, kepmask-ul or stacked-methanol)')

#area of circle spectra were extracted from
omega=pi*(0.4/3600.0*pi/180.0)**2

##### MEASURE LINE FLUXES
v_fitted=4.25


fitline='HDO_225'
restfreq_hdo=225.89672e9
restfreq_ch3ocho=225.90074e9
restfreq_13ch3cho=225.89294840e9
linefreqs=[restfreq_hdo,restfreq_ch3ocho,restfreq_13ch3cho]
lines=['           HDO','    CH$_3$OCHO','$^{13}$CH$_3$CHO']
lineamp=[0.14,0.1,0.075]
labels=['v','a1','a2','a3']
guesses=[v_fitted,3.0,2.4,0.5]
bounds=((v_fitted,v_fitted),(0.0,10.0),(0.0,10.0),(0.0,10.0))
rms_range=[225.89e9,225.9e9]
measure_range=[225.89e9,225.8975e9]
if measurement_type == 'H218O':
   # H218O model
   lineflux_225,e_lineflux_225=fit_HDO(fitline,'HDO_225_text_spectra.txt','H218O_text_spectral_template.txt',
               lines,linefreqs,lineamp,guesses,labels,bounds,rms_range,measure_range, template_file_coms='',useMCMC=False)
elif measurement_type == 'methanol':
   # CH3OH model
   lineflux_225,e_lineflux_225=fit_HDO(fitline,'HDO_225_text_spectra.txt','CH3OH_text_spectral_template.txt',
               lines,linefreqs,lineamp,guesses,labels,bounds,rms_range,measure_range, template_file_coms='',useMCMC=False)
elif measurement_type == 'mixed':
   # mixed model
   lineflux_225,e_lineflux_225=fit_HDO(fitline,'HDO_225_text_spectra.txt','H218O_text_spectral_template.txt',
               lines,linefreqs,lineamp,guesses,labels,bounds,rms_range,measure_range, template_file_coms='CH3OH_text_spectral_template.txt',useMCMC=False)
elif  measurement_type == 'stacked-methanol':
   lineflux_225,e_lineflux_225=fit_HDO(fitline,'HDO_225_stacked_text_spectra.txt','CH3OH_stacked_text_spectral_template.txt',
               lines,linefreqs,lineamp,guesses,labels,bounds,rms_range,measure_range, template_file_coms='CH3OH_stacked_text_spectral_template.txt',useMCMC=False)


fitline='HDO_241'
restfreq_hdo=241.56155e9
restfreq_ch3cho=241.56322770e9
restfreq_ch3cooh=241.56383270e9
restfreq_ch32co=241.55988690e9
linefreqs=[restfreq_hdo,restfreq_ch3cho,restfreq_ch32co]
lines=['             HDO','                CH$_3$CHO','(CH$_3$)$_2$CO']
lineamp=[0.19,0.19,0.19]
labels=['v','a1','a2','a3']
guesses=[v_fitted,3.0,0.5,1.5]
bounds=((v_fitted,v_fitted),(0.0,10.0),(0.0,10.0),(0.0,10.0))


rms_range=[241.555e9,241.564e9]
measure_range=[241.555e9,241.562e9]
if measurement_type == 'H218O':
   # H218O model
   lineflux_241,e_lineflux_241=fit_HDO(fitline,'HDO_241_text_spectra.txt','H218O_text_spectral_template.txt',
               lines,linefreqs,lineamp,guesses,labels,bounds,rms_range,measure_range, template_file_coms='',useMCMC=False)
elif measurement_type == 'methanol':
   # CH3OH model
   lineflux_241,e_lineflux_241=fit_HDO(fitline,'HDO_241_text_spectra.txt','CH3OH_text_spectral_template.txt',
               lines,linefreqs,lineamp,guesses,labels,bounds,rms_range,measure_range, template_file_coms='',useMCMC=False)
elif measurement_type == 'mixed':
   # mixed model
   lineflux_241,e_lineflux_241=fit_HDO(fitline,'HDO_241_text_spectra.txt','H218O_text_spectral_template.txt',
               lines,linefreqs,lineamp,guesses,labels,bounds,rms_range,measure_range, template_file_coms='CH3OH_text_spectral_template.txt',useMCMC=False)
elif  measurement_type == 'stacked-methanol':
   lineflux_241,e_lineflux_241=fit_HDO(fitline,'HDO_241_stacked_text_spectra.txt','CH3OH_stacked_text_spectral_template.txt',
               lines,linefreqs,lineamp,guesses,labels,bounds,rms_range,measure_range, template_file_coms='CH3OH_stacked_text_spectral_template.txt',useMCMC=False)

fitline='H218O_203'
restfreq_h218o=203.40752000e9
restfreq_ch3och3_1=203.41011260e9
restfreq_ch3och3_2=203.41143090e9
restfreq_ch3och3_3=203.40277790e9
linefreqs=[restfreq_h218o,restfreq_ch3och3_1,restfreq_ch3och3_2,restfreq_ch3och3_3]
lines=['H$_2$$^{18}$O','CH$_3$OCH$_3$','','CH$_3$OCH$_3$']
lineamp=[0.06,0.121,0.121,0.07]
labels=['v','a1','a2','a3','a4']
guesses=[v_fitted,0.1,0.1,0.1,0.1]
bounds=((v_fitted,v_fitted),(0.0,10.0),(0.0,10.0),(0.0,10.0),(0.0,10.0))
rms_range=[203.401e9,203.408e9]
measure_range=[203.401e9,203.408e9]


if measurement_type == 'H218O':
   # H218O model
   lineflux_203,e_lineflux_203=fit_HDO(fitline,'H218O_203_text_spectra.txt','H218O_text_spectral_template.txt',
               lines,linefreqs,lineamp,guesses,labels,bounds,rms_range,measure_range, template_file_coms='',useMCMC=False,label_v_offset=1.5,label_y_offset=0.01)
elif measurement_type == 'methanol':
   # CH3OH model
   lineflux_203,e_lineflux_203=fit_HDO(fitline,'H218O_203_text_spectra.txt','CH3OH_text_spectral_template.txt',
               lines,linefreqs,lineamp,guesses,labels,bounds,rms_range,measure_range, template_file_coms='',useMCMC=False)
elif measurement_type == 'mixed':
   # mixed model
   lineflux_203,e_lineflux_203=fit_HDO(fitline,'H218O_203_text_spectra.txt','H218O_text_spectral_template.txt',
               lines,linefreqs,lineamp,guesses,labels,bounds,rms_range,measure_range, template_file_coms='CH3OH_text_spectral_template.txt',useMCMC=False)
elif  measurement_type == 'stacked-methanol':
   # mixed model
   lineflux_203,e_lineflux_203=fit_HDO(fitline,'H218O_203_stacked_text_spectra.txt','CH3OH_stacked_text_spectral_template.txt',
               lines,linefreqs,lineamp,guesses,labels,bounds,rms_range,measure_range, template_file_coms='CH3OH_stacked_text_spectral_template.txt',useMCMC=False)



if measurement_type == 'kepmask':
   flux_dict=pickle.load(open("flux_dict.pickle", "rb"))
   # Keplerian mask measurements
   lineflux_225,e_lineflux_225=flux_dict['HDO_225']['flux_kepmask_raw'],flux_dict['HDO_225']['e_flux_kepmask_raw']
   #Keplerian mask measurements with estimated corrections for spectral range line flux is summed over
   lineflux_241,e_lineflux_241=flux_dict['HDO_241']['flux_kepmask_raw']+0.02*flux_dict['HDO_241']['flux_kepmask_raw'],flux_dict['HDO_241']['e_flux_kepmask_raw']
   lineflux_203,e_lineflux_203=flux_dict['H218O_203']['flux_kepmask_raw']+0.02*flux_dict['H218O_203']['flux_kepmask_raw'],flux_dict['H218O_203']['e_flux_kepmask_raw']


if measurement_type == 'kepmask-5-to-7':
   flux_dict=pickle.load(open("flux_dict.pickle", "rb"))
   # Keplerian mask measurements
   lineflux_225,e_lineflux_225=flux_dict['HDO_225']['flux_kepmask_5_to_7'],flux_dict['HDO_225']['e_flux_kepmask_5_to_7']
   lineflux_241,e_lineflux_241=flux_dict['HDO_241']['flux_kepmask_5_to_7'],flux_dict['HDO_241']['e_flux_kepmask_5_to_7']
   lineflux_203,e_lineflux_203=flux_dict['H218O_203']['flux_kepmask_5_to_7'],flux_dict['H218O_203']['e_flux_kepmask_5_to_7']

if measurement_type == 'kepmask-5-to-7-scale':
   flux_dict=pickle.load(open("flux_dict.pickle", "rb"))
   # Keplerian mask measurements
   lineflux_225,e_lineflux_225=flux_dict['HDO_225']['flux_kepmask_raw'],flux_dict['HDO_225']['e_flux_kepmask_raw']
   lineflux_241,e_lineflux_241=flux_dict['HDO_241']['flux_kepmask_scale_5_to_7_minus'],flux_dict['HDO_241']['e_flux_kepmask_scale_5_to_7_minus']
   lineflux_203,e_lineflux_203=flux_dict['H218O_203']['flux_kepmask_scale_5_to_7'],flux_dict['H218O_203']['e_flux_kepmask_scale_5_to_7']


if measurement_type == 'kepmask-ul':
   flux_dict=pickle.load(open("flux_dict.pickle", "rb"))
   lineflux_225,e_lineflux_225=flux_dict['HDO_225']['flux_kepmask_raw'],flux_dict['HDO_225']['e_flux_kepmask_raw']
   lineflux_241,e_lineflux_241=flux_dict['HDO_241']['flux_kepmask_scale_5_to_7_minus'],flux_dict['HDO_241']['e_flux_kepmask_scale_5_to_7_minus']
   lineflux_203,e_lineflux_203=flux_dict['H218O_203']['e_flux_kepmask_scale_5_to_7'],flux_dict['H218O_203']['e_flux_kepmask_scale_5_to_7']

if measurement_type == 'mom0':
   flux_dict=pickle.load(open("flux_dict.pickle", "rb"))
   #
   lineflux_225,e_lineflux_225=flux_dict['HDO_225']['flux']/dv_225,flux_dict['HDO_225']['e_flux']/dv_225
   lineflux_241,e_lineflux_241=flux_dict['HDO_241']['flux']/dv_241,flux_dict['HDO_241']['e_flux']/dv_241
   lineflux_203,e_lineflux_203=flux_dict['H218O_203']['flux']/dv_203,flux_dict['H218O_203']['e_flux']/dv_203
   #lineflux_203,e_lineflux_203=(0.09+0.1*0.09)/dv_203,0.01/dv_203

if measurement_type == 'mom0-scale':
   flux_dict=pickle.load(open("flux_dict.pickle", "rb"))
   lineflux_225,e_lineflux_225=flux_dict['HDO_225']['flux']/flux_dict['HDO_225']['flux_frac_5_to_7']/dv_225,flux_dict['HDO_225']['e_flux']/flux_dict['HDO_225']['flux_frac_5_to_7']/dv_225
   lineflux_241,e_lineflux_241=flux_dict['HDO_241']['flux']/flux_dict['HDO_225']['flux_frac_5_to_7']/dv_241,flux_dict['HDO_241']['e_flux']/flux_dict['HDO_225']['flux_frac_5_to_7']/dv_241
   lineflux_203,e_lineflux_203=flux_dict['H218O_203']['flux']/flux_dict['HDO_225']['flux_frac_5_to_7']/dv_203,flux_dict['H218O_203']['e_flux']/flux_dict['HDO_225']['flux_frac_5_to_7']/dv_203

fact=0.4
print(lineflux_225*fact,lineflux_241*fact,lineflux_203*fact)
print(e_lineflux_225*fact,e_lineflux_241*fact,e_lineflux_203*fact)
##### END MEASURE LINE FLUXES

plt.rcParams.update({'font.size': 10})

#calculate Nu
N_u_241_disk_avg_K=N_u_K(lineflux_241*1.0e5*dv_241*jy_to_k_aperture(0.4/3600.0,0.4/3600.0,nu_241),nu_241,A_ul_241)
N_u_225_disk_avg_K=N_u_K(lineflux_225*1.0e5*dv_225*jy_to_k_aperture(0.4/3600.0,0.4/3600.0,nu_225),nu_225,A_ul_225)
N_u_203_disk_avg_K=N_u_K(lineflux_203*1.0e5*dv_203*jy_to_k_aperture(0.4/3600.0,0.4/3600.0,nu_203),nu_203,A_ul_203)

# Nu uncertainties
e_N_u_241_disk_avg_K=N_u_241_disk_avg_K*e_lineflux_225/lineflux_225
e_N_u_225_disk_avg_K=N_u_225_disk_avg_K*e_lineflux_241/lineflux_241
e_N_u_203_disk_avg_K=N_u_203_disk_avg_K*e_lineflux_203/lineflux_203


# calculate line ratio and estimate Temperature
line_ratio_disk_avg=N_u_225_disk_avg_K*g_u_241/(N_u_241_disk_avg_K*g_u_225)
T_rot_diskavg=-(E_u_225-E_u_241)/(np.log(line_ratio_disk_avg))

#calculate uncertainty in Trot
e_N_u_225_d_N_u_241=N_u_225_disk_avg_K/(N_u_241_disk_avg_K)*((e_N_u_225_disk_avg_K/N_u_225_disk_avg_K)**2+(e_N_u_241_disk_avg_K/N_u_241_disk_avg_K)**2)**0.5
e_log_N_u_225_d_N_u_241=e_N_u_225_d_N_u_241/(N_u_225_disk_avg_K/(N_u_241_disk_avg_K))
e_Trot_diskavg=T_rot_diskavg*((e_log_N_u_225_d_N_u_241/np.log(line_ratio_disk_avg))**2)**0.5

Trot_disk_avg_string='{:0.0f}'.format(T_rot_diskavg)
# calculate Ntot given Trot and calculation of Qrot
N_241_disk_avg_K={'75 K': {'temp': 75.0, 'value': 0.0},
                  '180 K': {'temp': 180.0, 'value': 0.0},
                  'calcTrot': {'temp': T_rot_diskavg, 'value': 0.0},
                  'Tplus': {'temp': T_rot_diskavg+e_Trot_diskavg, 'value': 0.0},
                  'Tminus': {'temp': T_rot_diskavg-e_Trot_diskavg, 'value': 0.0},
                  '225 K': {'temp': 225.0, 'value': 0.0}}

N_225_disk_avg_K=copy.deepcopy(N_241_disk_avg_K)
N_203_disk_avg_K=copy.deepcopy(N_241_disk_avg_K) 
N_H2O_disk_avg_K=copy.deepcopy(N_241_disk_avg_K) 
N_HDO_disk_avg_K=copy.deepcopy(N_241_disk_avg_K) 
e_N_H2O_disk_avg_K=copy.deepcopy(N_H2O_disk_avg_K)
e_N_HDO_disk_avg_K=copy.deepcopy(N_H2O_disk_avg_K)
e_N_203_disk_avg_K=copy.deepcopy(N_203_disk_avg_K)
e_N_225_disk_avg_K=copy.deepcopy(N_225_disk_avg_K)
e_N_241_disk_avg_K=copy.deepcopy(N_241_disk_avg_K)
HDO_H2O_ratio_disk_avg=copy.deepcopy(N_241_disk_avg_K)
e_HDO_H2O_ratio_disk_avg=copy.deepcopy(N_241_disk_avg_K)

for key in N_241_disk_avg_K.keys():
   N_241_disk_avg_K[key]['value']=N_u_241_disk_avg_K/g_u_241*Qrot_HDO(N_241_disk_avg_K[key]['temp'])*np.exp(E_u_241/N_241_disk_avg_K[key]['temp'])
   N_225_disk_avg_K[key]['value']=N_u_225_disk_avg_K/g_u_225*Qrot_HDO(N_225_disk_avg_K[key]['temp'])*np.exp(E_u_225/N_225_disk_avg_K[key]['temp'])
   N_203_disk_avg_K[key]['value']=N_u_203_disk_avg_K/g_u_203*Qrot_H218O(N_203_disk_avg_K[key]['temp'])*np.exp(E_u_203/N_203_disk_avg_K[key]['temp'])
   N_H2O_disk_avg_K[key]['value']=N_u_203_disk_avg_K/g_u_203*Qrot_H218O(N_H2O_disk_avg_K[key]['temp'])*np.exp(E_u_203/N_203_disk_avg_K[key]['temp'])*560.0

   #calculate uncerts from line flux uncertainties
   e_N_241_disk_avg_K[key]['value']=N_241_disk_avg_K[key]['value']*((e_lineflux_241/lineflux_241)**2)**0.5
   e_N_225_disk_avg_K[key]['value']=N_225_disk_avg_K[key]['value']*((e_lineflux_225/lineflux_225)**2)**0.5
   e_N_203_disk_avg_K[key]['value']=N_203_disk_avg_K[key]['value']*((e_lineflux_203/lineflux_203)**2+(0.05)**2)**0.5
   #create error-weighted average N_HDO; formally the values will be the same since the should be identical
   if 'Tplus' not in key and 'Tminus' not in key:
      N_HDO_disk_avg_K[key]['value']=(N_241_disk_avg_K[key]['value']*1/e_N_241_disk_avg_K[key]['value']**2+N_225_disk_avg_K[key]['value']*1/e_N_225_disk_avg_K[key]['value']**2)/(1/e_N_241_disk_avg_K[key]['value']**2+1/e_N_225_disk_avg_K[key]['value']**2)
   else:
      N_HDO_disk_avg_K[key]['value']=(N_241_disk_avg_K[key]['value']*1/e_N_241_disk_avg_K['calcTrot']['value']**2+N_225_disk_avg_K[key]['value']*1/e_N_225_disk_avg_K['calcTrot']['value']**2)/(1/e_N_241_disk_avg_K['calcTrot']['value']**2+1/e_N_225_disk_avg_K['calcTrot']['value']**2)
      #N_HDO_disk_avg_K_Tplus=        (N_241_disk_avg_K_Tplus*       1/e_N_241_disk_avg_K**2+                      N_225_disk_avg_K_Tplus*        1/e_N_225_disk_avg_K**2)                     /(1/e_N_241_disk_avg_K**2                     +1/e_N_225_disk_avg_K**2)
   if key != 'calcTrot':
      e_N_H2O_disk_avg_K[key]['value']=N_H2O_disk_avg_K[key]['value']*((25.0/560.0)**2+(e_lineflux_203/lineflux_203)**2+(0.05)**2)**0.5
      e_N_HDO_disk_avg_K[key]['value']=((e_N_241_disk_avg_K[key]['value']**2+e_N_225_disk_avg_K[key]['value']**2)**0.5/2.0+ (N_HDO_disk_avg_K[key]['value']*0.05)**2)**0.5 

#calculate uncertainty contribution of Trot uncertainty to H2O column using 1-sigma error of T_rot
N_H2O_disk_avg_Terror_comp=(abs(N_H2O_disk_avg_K['Tplus']['value']-N_H2O_disk_avg_K['calcTrot']['value'])/N_H2O_disk_avg_K['calcTrot']['value']+abs(N_H2O_disk_avg_K['calcTrot']['value']-N_H2O_disk_avg_K['Tminus']['value'])/N_H2O_disk_avg_K['calcTrot']['value'])/2.0
N_HDO_disk_avg_Terror_comp=(abs(N_HDO_disk_avg_K['Tplus']['value']-N_HDO_disk_avg_K['calcTrot']['value'])/N_HDO_disk_avg_K['calcTrot']['value']+abs(N_HDO_disk_avg_K['calcTrot']['value']-N_HDO_disk_avg_K['Tminus']['value'])/N_HDO_disk_avg_K['calcTrot']['value'])/2.0

#calculate uncertainty of H2O due to temperature error using 1-sigma error of T_rot
e_N_H2O_disk_avg_K['calcTrot']['value']=N_H2O_disk_avg_K['calcTrot']['value']*((25.0/560.0)**2+(e_lineflux_203/lineflux_203)**2+(0.05)**2+(N_H2O_disk_avg_Terror_comp)**2)**0.5
e_N_HDO_disk_avg_K['calcTrot']['value']=N_HDO_disk_avg_K['calcTrot']['value']*((25.0/560.0)**2+(e_lineflux_203/lineflux_203)**2+(0.05)**2+(N_H2O_disk_avg_Terror_comp)**2)**0.5
e_N_HDO_disk_avg_K['calcTrot']['value']=((e_N_241_disk_avg_K['calcTrot']['value']**2+e_N_225_disk_avg_K['calcTrot']['value']**2)**0.5/2.0+ (N_HDO_disk_avg_K['calcTrot']['value']*0.05)**2 + (N_HDO_disk_avg_K['calcTrot']['value']*N_HDO_disk_avg_Terror_comp)**2)**0.5 # with flux denisty uncertainty

for key in HDO_H2O_ratio_disk_avg.keys():
   HDO_H2O_ratio_disk_avg[key]['value']=N_HDO_disk_avg_K[key]['value']/N_H2O_disk_avg_K[key]['value']
   e_HDO_H2O_ratio_disk_avg[key]['value']=HDO_H2O_ratio_disk_avg[key]['value']*((e_N_HDO_disk_avg_K[key]['value']/N_HDO_disk_avg_K[key]['value'])**2+(e_N_H2O_disk_avg_K[key]['value']/N_H2O_disk_avg_K[key]['value'])**2)**0.5

print('Disk Avg. T_rot: {:0.1f}+/-{:0.1f}'.format(T_rot_diskavg,e_Trot_diskavg))
print('Q_HDO({:0.1f}) = {:0.1f} '.format(T_rot_diskavg,Qrot_HDO(T_rot_diskavg)))
print('Q_H218O({:0.1f}) = {:0.1f} '.format(T_rot_diskavg,Qrot_H218O(T_rot_diskavg)))
print('Disk Avg N_tot HDO 225 GHz: {:0.3e}+/-{:0.3e}'.format(N_225_disk_avg_K['calcTrot']['value'],e_N_225_disk_avg_K['calcTrot']['value']))
print('Disk Avg N_tot HDO 241 GHz:  {:0.3e}+/-{:0.3e}'.format(N_241_disk_avg_K['calcTrot']['value'],e_N_241_disk_avg_K['calcTrot']['value']))
print('Disk Avg N_tot HDO Avg:  {:0.3e}+/-{:0.3e}'.format(N_HDO_disk_avg_K['calcTrot']['value'],e_N_HDO_disk_avg_K['calcTrot']['value']))
print('Disk Avg N_tot H218O 203 GHz: {:0.3e}+/-{:0.3e}'.format(N_203_disk_avg_K['calcTrot']['value'],e_N_203_disk_avg_K['calcTrot']['value']))
print('Inferred Disk Avg N_tot H2O: {:0.3e}+/-{:0.3e}'.format(N_H2O_disk_avg_K['calcTrot']['value'],e_N_H2O_disk_avg_K['calcTrot']['value']))
print('HDO/H2O ratio (T={:0.1f}K): {:0.3e}+/-{:0.3e}'.format(T_rot_diskavg,HDO_H2O_ratio_disk_avg['calcTrot']['value'],e_HDO_H2O_ratio_disk_avg['calcTrot']['value']))
print('Water mass: {:0.3e} g'.format(N_H2O_disk_avg_K['calcTrot']['value']*3.14*160**2*1.496e13**2*1.67e-24*(1.0+1.0+16.0)))
print('Earth Oceans: {:0.3e} '.format(N_H2O_disk_avg_K['calcTrot']['value']*3.14*160**2*1.496e13**2*1.67e-24*(1.0+1.0+16.0)/1.384e24))
for key in N_241_disk_avg_K.keys():
   if 'Tplus' not in key and 'Tminus' not in key:
      print('HDO/H2O ratio (T={}): {:0.3e}+/-{:0.3e}'.format(key, HDO_H2O_ratio_disk_avg[key]['value'],e_HDO_H2O_ratio_disk_avg[key]['value']))


##### BEGIN RADIAL PROFILE ANALYSIS HERE

#ignore zero/inf divide warnings
np.seterr(all='ignore')

#radial_profiles=pickle.load(open("radial_profiles_Jykms.pickle", "rb"))
radial_profiles_Kkms=pickle.load(open("radial_profiles_Kkms.pickle", "rb"))
#for i in range(len(radial_profiles['HDO 225.535GHz']['x'])):
#   print(radial_profiles['HDO 225.535GHz']['x'][i],radial_profiles['HDO 241.558 GHz']['x'][i])

#spatial coordinates are slightly different # so need to interpolate
#HDO_241_profile='HDO 241.558 GHz'
#HDO_225_profile='HDO 225.535GHz'
# use profiles from data convolved to same beam as H218O
HDO_225_profile='HDO 225.535GHz'
HDO_241_profile='HDO 241.558 GHz'

hdo_225_radial_profile_interp_Kkms=np.interp(radial_profiles_Kkms[HDO_241_profile]['x'],radial_profiles_Kkms[HDO_225_profile]['x'],radial_profiles_Kkms[HDO_225_profile]['y'])
e_hdo_225_radial_profile_interp_Kkms=np.interp(radial_profiles_Kkms[HDO_241_profile]['x'],radial_profiles_Kkms[HDO_225_profile]['x'],radial_profiles_Kkms[HDO_225_profile]['dy'])
h218o_203_radial_profile_interp_Kkms=np.interp(radial_profiles_Kkms[HDO_241_profile]['x'],radial_profiles_Kkms['H218O_203']['x'],radial_profiles_Kkms['H218O_203']['y'])
e_h218o_203_radial_profile_interp_Kkms=np.interp(radial_profiles_Kkms[HDO_241_profile]['x'],radial_profiles_Kkms['H218O_203']['x'],radial_profiles_Kkms['H218O_203']['dy'])


N_u_241_K=N_u_K(radial_profiles_Kkms[HDO_241_profile]['y']*1.0e5,nu_241,A_ul_241)
N_u_225_K=N_u_K(hdo_225_radial_profile_interp_Kkms*1.0e5,nu_225,A_ul_225)
N_u_203_K=N_u_K(h218o_203_radial_profile_interp_Kkms*1.0e5,nu_203,A_ul_203)
e_N_u_241_K=N_u_241_K*radial_profiles_Kkms[HDO_241_profile]['dy']/radial_profiles_Kkms[HDO_241_profile]['y']
e_N_u_225_K=N_u_225_K*e_hdo_225_radial_profile_interp_Kkms/hdo_225_radial_profile_interp_Kkms
e_N_u_203_K=N_u_203_K*e_h218o_203_radial_profile_interp_Kkms/h218o_203_radial_profile_interp_Kkms
line_ratio_profile=(N_u_225_K*g_u_241)/(N_u_241_K*g_u_225)
T_rot_profile=-(E_u_225-E_u_241)/(np.log(line_ratio_profile))


e_N_u_225_d_N_u_241=N_u_225_K/(N_u_241_K)*((e_N_u_225_K/N_u_225_K)**2+(e_N_u_241_K/N_u_241_K)**2)**0.5
e_log_N_u_225_d_N_u_241=e_N_u_225_d_N_u_241/(N_u_225_K/(N_u_241_K))
e_T_rot_profile=T_rot_profile*((e_log_N_u_225_d_N_u_241/np.log(line_ratio_profile))**2)**0.5

N_241_disk_profile={'calcTrot': {'temp': T_rot_diskavg},
                  '75 K': {'temp': 75.0},
                  '225 K': {'temp': 225.0},
                  'calcTrot_profile': {'temp': T_rot_profile},
                  'Tplus_profile': {'temp': T_rot_profile+e_T_rot_profile},
                  'Tminus_profile': {'temp': T_rot_profile-e_T_rot_profile},
                  }

N_225_disk_profile=copy.deepcopy(N_241_disk_profile)
N_203_disk_profile=copy.deepcopy(N_241_disk_profile) 
N_H2O_disk_profile=copy.deepcopy(N_241_disk_profile) 
N_HDO_disk_profile=copy.deepcopy(N_241_disk_profile) 
e_N_H2O_disk_profile=copy.deepcopy(N_H2O_disk_profile)
e_N_HDO_disk_profile=copy.deepcopy(N_H2O_disk_profile)
e_N_203_disk_profile=copy.deepcopy(N_203_disk_profile)
e_N_225_disk_profile=copy.deepcopy(N_225_disk_profile)
e_N_241_disk_profile=copy.deepcopy(N_241_disk_profile)
HDO_H2O_ratio_disk_profile=copy.deepcopy(N_241_disk_profile)
e_HDO_H2O_ratio_disk_profile=copy.deepcopy(N_241_disk_profile)

for key in N_241_disk_profile.keys():
   N_241_disk_profile[key]['profile']=N_u_241_K/g_u_241*Qrot_HDO(N_241_disk_profile[key]['temp'])*np.exp(E_u_241/N_241_disk_profile[key]['temp'])
   N_225_disk_profile[key]['profile']=N_u_225_K/g_u_225*Qrot_HDO(N_241_disk_profile[key]['temp'])*np.exp(E_u_225/N_225_disk_profile[key]['temp'])
   N_203_disk_profile[key]['profile']=N_u_203_K/g_u_203*Qrot_H218O(N_203_disk_profile[key]['temp'])*np.exp(E_u_203/N_203_disk_profile[key]['temp'])
   N_H2O_disk_profile[key]['profile']=N_203_disk_profile[key]['profile']*560.0
   e_N_241_disk_profile[key]['profile']=N_241_disk_profile[key]['profile']*radial_profiles_Kkms[HDO_241_profile]['dy']/radial_profiles_Kkms[HDO_241_profile]['y']
   e_N_225_disk_profile[key]['profile']=N_225_disk_profile[key]['profile']*e_hdo_225_radial_profile_interp_Kkms/hdo_225_radial_profile_interp_Kkms
   e_N_203_disk_profile[key]['profile']=N_203_disk_profile[key]['profile']*e_h218o_203_radial_profile_interp_Kkms/h218o_203_radial_profile_interp_Kkms
   e_N_H2O_disk_profile[key]['profile']=N_203_disk_profile[key]['profile']*e_h218o_203_radial_profile_interp_Kkms/h218o_203_radial_profile_interp_Kkms*560.0
   N_HDO_disk_profile[key]['profile']=(N_241_disk_profile[key]['profile']*1/e_N_241_disk_profile[key]['profile']**2+N_225_disk_profile[key]['profile']*1/e_N_225_disk_profile[key]['profile']**2)/(1/e_N_241_disk_profile[key]['profile']**2+1/e_N_225_disk_profile[key]['profile']**2)
   #only use HDO 225 GHz since 241 is quite contaminated in this implementation
   #N_HDO_disk_profile[key]['profile']=N_225_disk_profile[key]['profile']
   if key == 'calcTrot':
      e_N_HDO_disk_profile[key]['profile']=((2.0/(1/e_N_241_disk_profile[key]['profile']**2+1/e_N_225_disk_profile[key]['profile']**2))+ (N_HDO_disk_profile[key]['profile']*0.05)**2)**0.5
      #only HDO 225 again      
      #e_N_HDO_disk_profile[key]['profile']=((e_N_225_disk_profile[key]['profile']**2)+ (N_HDO_disk_profile[key]['profile']*0.05)**2)**0.5
      e_N_H2O_disk_profile[key]['profile']=N_H2O_disk_profile[key]['profile']*((25.0/560.0)**2+(e_N_203_disk_profile[key]['profile']/N_203_disk_profile[key]['profile'])**2+(0.05)**2+(0.044)**2)**0.5

N_HDO_radial_profile_avg_Trot_profile_Terror_comp=((N_HDO_disk_profile['Tplus_profile']['profile']-N_HDO_disk_profile['calcTrot_profile']['profile'])/N_HDO_disk_profile['calcTrot_profile']['profile']+(N_HDO_disk_profile['Tminus_profile']['profile']-N_HDO_disk_profile['calcTrot_profile']['profile'])/N_HDO_disk_profile['calcTrot_profile']['profile'])/2.0
e_N_HDO_disk_profile['calcTrot_profile']['profile']=((2.0/(1/e_N_241_disk_profile['calcTrot_profile']['profile']**2+1/e_N_225_disk_profile['calcTrot_profile']['profile']**2))+ (N_HDO_disk_profile['calcTrot_profile']['profile']*0.05)**2 + (N_HDO_disk_profile['calcTrot_profile']['profile']*N_HDO_radial_profile_avg_Trot_profile_Terror_comp)**2)**0.5

N_H2O_radial_profile_avg_Trot_profile_Terror_comp=((N_H2O_disk_profile['Tplus_profile']['profile']-N_H2O_disk_profile['calcTrot_profile']['profile'])/N_H2O_disk_profile['calcTrot_profile']['profile']+(N_H2O_disk_profile['Tminus_profile']['profile']-N_H2O_disk_profile['calcTrot_profile']['profile'])/N_H2O_disk_profile['calcTrot_profile']['profile'])/2.0
e_N_H2O_disk_profile['calcTrot_profile']['profile']=N_H2O_disk_profile['calcTrot_profile']['profile']*((25.0/560.0)**2+(e_N_203_disk_profile['calcTrot_profile']['profile']/N_203_disk_profile['calcTrot_profile']['profile'])**2+(0.05)**2 + N_H2O_radial_profile_avg_Trot_profile_Terror_comp**2)**0.5

for key in HDO_H2O_ratio_disk_profile.keys():
   HDO_H2O_ratio_disk_profile[key]['profile']=N_HDO_disk_profile[key]['profile']/N_H2O_disk_profile[key]['profile']
   if ('profile' in e_N_HDO_disk_profile[key].keys()) and ('profile' in e_N_H2O_disk_profile[key].keys()):
      e_HDO_H2O_ratio_disk_profile[key]['profile']=HDO_H2O_ratio_disk_profile[key]['profile']*((e_N_HDO_disk_profile[key]['profile']/N_HDO_disk_profile[key]['profile'])**2+(e_N_H2O_disk_profile[key]['profile']/N_H2O_disk_profile[key]['profile'])**2)**0.5



fig, ax = plt.subplots()
maxy=np.nanmax(T_rot_diskavg)
print('Max Trot: ',maxy)
ax.errorbar(radial_profiles_Kkms[HDO_241_profile]['x'], T_rot_profile, T_rot_profile/100.0, capsize=1.25, capthick=1.25, lw=1.25,label='T$_{rot}$ HDO')
ax.legend()
ax.set_xlim(np.min(radial_profiles_Kkms[HDO_241_profile]['x']), 0.45)
ax.set_xlabel('Radius (arcsec)')
ax.set_ylabel('T$_{rot}$ (K)')
ax.set_title('HDO T$_{rot}$ Profile')
plt.savefig('HDO_T_rot_Radial_Profile.png',dpi=200.0)
plt.savefig('HDO_T_rot_Radial_Profile.pdf',dpi=200.0)


line_ratio_radial_profile=radial_profiles_Kkms[HDO_241_profile]['y']/hdo_225_radial_profile_interp_Kkms
e_line_ratio_radial_profile=line_ratio_radial_profile*((radial_profiles_Kkms[HDO_241_profile]['dy']/radial_profiles_Kkms[HDO_241_profile]['y'])**2+(e_hdo_225_radial_profile_interp_Kkms/hdo_225_radial_profile_interp_Kkms)**2)**0.5


fig, ax = plt.subplots()
maxy=np.nanmax(T_rot_diskavg)
ax.errorbar(radial_profiles_Kkms[HDO_241_profile]['x'], line_ratio_radial_profile,e_line_ratio_radial_profile, capsize=1.25, capthick=1.25, lw=1.25,label='T$_{rot}$ HDO')
ax.legend()
ax.set_xlim(np.min(radial_profiles_Kkms[HDO_241_profile]['x']), 0.45)
ax.set_xlabel('Radius (arcsec)')
ax.set_ylabel('T$_{rot}$ (K)')
ax.set_title('HDO Line Ratio Profile')
plt.savefig('HDO_Line_Ratios.png',dpi=200.0)
plt.savefig('HDO_Line_Ratios.pdf',dpi=200.0)


fig, ax = plt.subplots()
ax.errorbar(radial_profiles_Kkms[HDO_241_profile]['x'][0:17], N_H2O_disk_profile['calcTrot_profile']['profile'][0:17]/560.,lw=1.25,label='H$_2$$^{18}$O',color='blue')
ax.errorbar(radial_profiles_Kkms[HDO_241_profile]['x'][0:17], N_H2O_disk_profile['calcTrot_profile']['profile'][0:17]/560., e_N_H2O_disk_profile['calcTrot_profile']['profile'][0:17]/560., capsize=1.25, capthick=1.25, lw=1.25 ,label='_H$_2$$^{18}$O 203.407 GHz',color='blue')
ax.errorbar(radial_profiles_Kkms[HDO_241_profile]['x'][0:17], N_HDO_disk_profile['calcTrot_profile']['profile'][0:17] ,lw=1.25,label='HDO',color='orange')
ax.errorbar(radial_profiles_Kkms[HDO_241_profile]['x'][0:17], N_HDO_disk_profile['calcTrot_profile']['profile'][0:17], e_N_HDO_disk_profile['calcTrot_profile']['profile'][0:17], capsize=1.25, capthick=1.25,lw=1.25,label='_HDO 225.896 and 241.561 GHz',color='orange')
ax.plot(radial_profiles_Kkms[HDO_241_profile]['x'][0:17],T_rot_profile[0:17],label='HDO T$_{ex}$ Profile',color='red', lw=1.25,linestyle='dashed')
ax.set_yscale('log')

ax.set_xlim(np.min(radial_profiles_Kkms[HDO_241_profile]['x']),0.45)
ax.set_xlabel('Radius (arcsec)',fontsize='large')
ax.set_ylabel('Column Density (cm$^{-2}$)',fontsize='large')
ax.set_ylim(0.5e14,1e17)
ax.set_xlim(0,0.425)
ax.tick_params(axis='x', labelsize='large')
ax.tick_params(axis='y', labelsize='large')
ax.fill_between([0,0.1], 1e12, 1e17,
                facecolor='gray', alpha=0.75)
ax.text(0.04,3e16,' Optically\n Thick\n Dust')
new_tick_locations = np.array([0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45])
ax2 = ax.twiny()
def tick_function(X):
    V = X*400.0
    return ["%.0f" % z for z in V]

ax2.set_xlim(ax.get_xlim())
ax2.set_xticks(new_tick_locations)
ax2.set_xticklabels(tick_function(new_tick_locations),fontsize='large')
ax2.set_xlabel("Radius (au)",fontsize='large')

new_tick_locations = np.array([1.0e14,1.0e15,1.0e16,1.0e17])
ax3 = ax.twinx()
def tick_function(X):
    #V = X*400.0
    V=X/1.0e14
    return ["%.0f" % z for z in V]

ax3.set_yscale('log')
ax3.set_ylim(ax.get_ylim())
#ax3.set_ylim([1.0,1000.0])
ax3.set_yticks(new_tick_locations)
ax3.set_yticklabels(tick_function(new_tick_locations),fontsize='large')
ax3.set_ylabel("Temperature(K)",fontsize='large',labelpad=-2)
ax3.errorbar(radial_profiles_Kkms[HDO_241_profile]['x'][0:17],T_rot_profile[0:17]*1.0e14,e_T_rot_profile[0:17]*1.0e14, capsize=1.25, capthick=1.25,label='_Temperature Profile',color='red',linestyle='dashed')

ax.errorbar([0.2,],[T_rot_diskavg*1.0e14],\
        [e_Trot_diskavg*1.0e14],\
        linestyle='none',marker='o',label='Disk Avg. Temp.',color='red')

ax.errorbar([0.215,],[N_HDO_disk_avg_K['calcTrot']['value']],\
        [e_N_HDO_disk_avg_K['calcTrot']['value']],\
        linestyle='none',marker='s',label='Disk Avg. HDO',color='orange')

ax.errorbar([0.2,],[N_203_disk_avg_K['calcTrot']['value']],\
        [e_N_203_disk_avg_K['calcTrot']['value']],\
        linestyle='none',marker='D',label='Disk Avg. H$_2$$^{18}$O',color='blue')

ax.text(0.4,6e16,'a)',fontsize='large')
ax.legend(loc='lower left',fontsize='large')

plt.savefig('HDO_H218O_Column_Density_Profile_K_Trot_profile.png',dpi=200.0)
plt.savefig('HDO_H218O_Column_Density_Profile_K_Trot_profile.pdf',dpi=200.0)



fig, ax = plt.subplots()
ax.errorbar(radial_profiles_Kkms[HDO_241_profile]['x'], N_241_disk_profile['calcTrot_profile']['profile'], e_N_241_disk_profile['calcTrot_profile']['profile'], capsize=1.25, capthick=1.25, lw=1.25,label='HDO 241.896 GHz')
ax.errorbar(radial_profiles_Kkms[HDO_241_profile]['x'], N_225_disk_profile['calcTrot_profile']['profile'], e_N_225_disk_profile['calcTrot_profile']['profile'], capsize=1.25, capthick=1.25, lw=1.25,label='HDO 225.561 GHz')
ax.errorbar(radial_profiles_Kkms[HDO_241_profile]['x'], N_203_disk_profile['calcTrot_profile']['profile'], e_N_203_disk_profile['calcTrot_profile']['profile'], capsize=1.25, capthick=1.25, lw=1.25,label='H$_2$$^{18}$O 203.407 GHz')
#ax.errorbar(radial_profiles[HDO_241_profile]['x'], N_203_K*560*3.0, N_203_K*560*3.0*((e_N_203_K/N_203_K)**2+(25/560)**2)**0.5, capsize=1.25, capthick=1.25, lw=1.25,label='H$_2$O Inferred')
ax.set_yscale('log')
ax.legend()
ax.set_xlim(np.min(radial_profiles_Kkms[HDO_241_profile]['x']),0.45)
ax.set_xlabel('Radius (arcsec)')
ax.set_ylabel('Column Density (cm$^{-2}$)')
ax.set_title('Water Column Density Profiles')
plt.savefig('HDO_H218O_Column_Density_Profile_K.png',dpi=200.0)
plt.savefig('HDO_H218O_Column_Density_Profile_K.pdf',dpi=200.0)

fig, ax = plt.subplots()
ax.errorbar(radial_profiles_Kkms[HDO_241_profile]['x'], N_203_disk_profile['calcTrot_profile']['profile'], e_N_203_disk_profile['calcTrot_profile']['profile'], capsize=1.25, capthick=1.25, lw=1.25,label='H$_2$$^{18}$O 203.40725 GHz (K)')
ax.set_yscale('log')
ax.legend()
ax.set_xlim(np.min(radial_profiles_Kkms[HDO_241_profile]['x']),0.45)
ax.set_xlabel('Radius (arcsec)')
ax.set_ylabel('Column Density (cm$^{-2}$)')
ax.set_title('H$_2$$^{18}$O Column Density Profile')
plt.savefig('H218O_Column_Density_Profile.png',dpi=200.0)
plt.savefig('H218O_Column_Density_Profile.pdf',dpi=200.0)


Trot_string='{:0.1f} K'.format(T_rot_diskavg)
fig, ax = plt.subplots() 
ax.plot(radial_profiles_Kkms[HDO_241_profile]['x'][0:17], HDO_H2O_ratio_disk_profile['calcTrot']['profile'][0:17],  lw=1.25,label='HDO/H$_2$$^{18}$O T$_{rot}$ = '+Trot_string,color='black',linestyle='--')
ax.errorbar(radial_profiles_Kkms[HDO_241_profile]['x'][0:17], HDO_H2O_ratio_disk_profile['calcTrot_profile']['profile'][0:17], e_HDO_H2O_ratio_disk_profile['calcTrot_profile']['profile'][0:17], capsize=1.25, capthick=1.25, lw=1.25,label='HDO/H$_2$$^{18}$O T$_{rot}$ Profile',color='black')
ax.errorbar([0.2,],[HDO_H2O_ratio_disk_avg['calcTrot']['value']],\
        [e_HDO_H2O_ratio_disk_avg['calcTrot']['value']],\
        linestyle='none',marker='o',label='HDO/H$_2$$^{18}$O Disk Averaged')

ax.fill_between([0,0.1], 0, 1,
                facecolor='gray', alpha=0.75)
ax.text(0.02,0.003,' Optically\n Thick\n Dust')
ax.tick_params(axis='x', labelsize='large')
ax.tick_params(axis='y', labelsize='large')

ax.legend(fontsize='large')
ax.set_xlim(np.min(radial_profiles_Kkms[HDO_241_profile]['x']), 0.45)
ax.set_xlim(0.0, 0.425)
ax.set_ylim(0.0001,0.01)
ax.set_xlabel('Radius (arcsec)',fontsize='large')
ax.set_ylabel('HDO/H$_2$O',fontsize='large')
ax.set_yscale('log')

new_tick_locations = np.array([0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45])
ax2 = ax.twiny()
def tick_function(X):
    V = X*400.0
    return ["%.0f" % z for z in V]

ax2.set_xlim(ax.get_xlim())
ax2.set_xticks(new_tick_locations)
ax2.set_xticklabels(tick_function(new_tick_locations),fontsize='large')
ax2.set_xlabel("Radius (au)",fontsize='large')
ax.text(0.4,7.3e-3,'b)',fontsize='large')
plt.savefig('HDO_H2O_Column_Density_Profile.png',dpi=200.0)
plt.savefig('HDO_H2O_Column_Density_Profile.pdf',dpi=200.0)

fig, ax = plt.subplots()

ax.errorbar([E_u_241,E_u_225],np.log([N_u_241_disk_avg_K/g_u_241,N_u_225_disk_avg_K/g_u_225]),\
            np.log([N_u_241_disk_avg_K/g_u_241,N_u_225_disk_avg_K/g_u_225])*0.01,\
            capsize=1.25, capthick=1.25,marker='o', linestyle='none',label='HDO Transitions')
ax.legend()
ax.set_xlim(0,200)
ax.set_xlabel('E/k (K)')
ax.set_ylabel('ln (N$_u$/g$_u$)')
ax.set_title('HDO Rotation Diagram')
plt.savefig('HDO_Rotation_diagram.png',dpi=200.0)
plt.savefig('HDO_Rotation_diagram.pdf',dpi=200.0)


max_HDO_column_density=np.nanmax(N_HDO_disk_profile['calcTrot_profile']['profile'][0:17])

slope=np.zeros(17)
radius=np.zeros(17)
for i in range(17):
   slope[i]=(N_HDO_disk_profile['calcTrot_profile']['profile'][i+1]-N_HDO_disk_profile['calcTrot_profile']['profile'][i])/(radial_profiles_Kkms[HDO_241_profile]['x'][i+1]-radial_profiles_Kkms[HDO_241_profile]['x'][i])
   radius[i]=radial_profiles_Kkms[HDO_241_profile]['x'][i]+(radial_profiles_Kkms[HDO_241_profile]['x'][i+1]-radial_profiles_Kkms[HDO_241_profile]['x'][i])/2
print(slope)
print(radius*d)


print(max_HDO_column_density)
index=np.where(N_HDO_disk_profile['calcTrot_profile']['profile'][0:17]==max_HDO_column_density)
print('Peak HDO column density Radius: {:0.2f} au'.format(radial_profiles_Kkms[HDO_241_profile]['x'][index[0][0]]*d))
index=np.where(N_HDO_disk_profile['calcTrot_profile']['profile'][0:17]<max_HDO_column_density*0.5)
print('50% HDO column density Radius: {:0.2f} au'.format(radial_profiles_Kkms[HDO_241_profile]['x'][index[0][0]]*d))


fig, ax = plt.subplots()
ax.errorbar(radial_profiles_Kkms[HDO_241_profile]['x'], N_241_disk_profile['calcTrot_profile']['profile'], e_N_241_disk_profile['calcTrot_profile']['profile'], capsize=1.25, capthick=1.25, lw=1.25,label='HDO 241.896 GHz')
ax.errorbar(radial_profiles_Kkms[HDO_241_profile]['x'], N_225_disk_profile['calcTrot_profile']['profile'], e_N_225_disk_profile['calcTrot_profile']['profile'], capsize=1.25, capthick=1.25, lw=1.25,label='HDO 225.561 GHz')
ax.errorbar(radial_profiles_Kkms[HDO_241_profile]['x'], N_203_disk_profile['calcTrot_profile']['profile'], e_N_203_disk_profile['calcTrot_profile']['profile'], capsize=1.25, capthick=1.25, lw=1.25,label='H$_2$$^{18}$O 203.407 GHz')
ax.plot(radius,np.abs(slope),  lw=1.25,label='HDO/H$_2$$^{18}$O T$_{rot}$ = '+Trot_string,color='cyan',linestyle='dotted')
#ax.errorbar(radial_profiles[HDO_241_profile]['x'], N_203_K*560*3.0, N_203_K*560*3.0*((e_N_203_K/N_203_K)**2+(25/560)**2)**0.5, capsize=1.25, capthick=1.25, lw=1.25,label='H$_2$O Inferred')
ax.set_yscale('log')
ax.legend()
ax.set_xlim(np.min(radial_profiles_Kkms[HDO_241_profile]['x']),0.45)
ax.set_xlabel('Radius (arcsec)')
ax.set_ylabel('Column Density (cm$^{-2}$)')
ax.set_title('Water Column Density Profiles')
plt.savefig('HDO_H218O_Column_Density_Profile_K_wslope.png',dpi=200.0)
plt.savefig('HDO_H218O_Column_Density_Profile_K_wslope.pdf',dpi=200.0)


